Contents

Advanced Data Visualization

In the last module we started building some visualizations in order to answer specific substantive questions. In this module we will look at generating more advanced graphics that have annotations, combine several graphics into one, build some interactive graphics, and also do some mapping.

Combining Plots with {patchwork}

Often you have to combine and place multiple graphics into a single canvas. There are a few ways to do this but the easiest way is that offered by the {patchwork} package. Let us use the diamonds data for this section, a data frame with 53940 rows and 10 variables:

Variable

Description

price

price in US dollars ($326–$18,823)

carat

weight of the diamond (0.2–5.01)

cut

quality of the cut (Fair, Good, Very Good, Premium, Ideal)

color

diamond colour, from D (best) to J (worst)

clarity

a measurement of how clear the diamond is (I1 (worst), SI2, SI1, VS2, VS1, VVS2, VVS1, IF (best))

x

length in mm (0–10.74)

y

width in mm (0–58.9)

z

depth in mm (0–31.8)

depth

total depth percentage = z / mean(x, y) = 2 * z / (x + y) (43–79)

table

width of top of diamond relative to widest point (43–95)

The Basics

To combine multiple plots, we need to save each plot with a unique name. I am calling them p1, p2, etc. Let us generate four plots, each different from all others.

#install.packages(c("patchwork", "tidylog", "leaflet", "htmltools", "highcharter"))

library(patchwork)
library(tidyverse)
library(tidylog)

── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──

✔ ggplot2 3.3.5     ✔ purrr   0.3.4
✔ tibble  3.1.6     ✔ dplyr   1.0.7
✔ tidyr   1.1.4     ✔ stringr 1.4.0
✔ readr   2.1.1     ✔ forcats 0.5.1

── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

Attaching package: ‘tidylog’

The following objects are masked from ‘package:dplyr’:

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup

The following objects are masked from ‘package:tidyr’:

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount

The following object is masked from ‘package:stats’:

    filter

options(repr.plot.width = 16, repr.plot.height = 8)

ggplot() +
  geom_bar(data = diamonds,
           aes(x = cut, fill = cut)
           ) +
  labs(x = "Cut of the Diamond", y = "Frequency") +
  theme(legend.position = "none") -> p1

ggplot() +
  geom_bar(data = diamonds,
           aes(x = color, fill = color)
           ) +
  labs(x = "Color of the Diamond", y = "Frequency") +
  theme(legend.position = "none") -> p2

ggplot() +
  geom_point(data = diamonds,
           aes(x = carat, y = price, color = cut)
           ) +
  labs(x = "Weight of the Diamond", y = "Price of the Diamond", color = "") +
  theme(legend.position = "bottom") -> p3

ggplot() +
  geom_boxplot(data = diamonds,
           aes(x = price, y = clarity, fill = cut)
           ) +
  labs(y = "Clarity of the Diamond", x = "Price of the Diamond", fill = "") +
  theme(legend.position = "bottom") -> p4

Let us see each plot in turn so we know what they look like.

p1; p2; p3; p4
../_images/Module06_10_0.png ../_images/Module06_10_1.png ../_images/Module06_10_2.png ../_images/Module06_10_3.png

Now we combine p1 through p3 on a single canvas

p1 + p2 + p3
../_images/Module06_12_0.png

Notice the default layout here: p1 + p2 + p3 gives us the plots all in a row.

But you may have other plans, for example, to put the scatterplot in a row all its own.

(p1 + p2) / p3
../_images/Module06_14_0.png

Now we have p3 in the second row, all by itself. Note that this was achieved via the / operator and by coercing p1 and p2 into a single row via (p1 + p2).

What if we used | instead?

p1 | (p2 + p3)
../_images/Module06_16_0.png

You ended up with two columns, the first containing only p1 and the second containing p2 and p3.

Note the difference between | and /. For example, note the following setup:

p1 | (p2 / p3)
../_images/Module06_18_0.png

You get p1 in one-half of the canvas, and then p2 and p3 are split into two rows in the remaining half of the canvas.

What if we wanted to squeeze in the fourth plot?

options(repr.plot.width = 16, repr.plot.height = 16)

(p1 + p2) / (p3 + p4)
../_images/Module06_20_0.png

Here you asked for p1 and p2 to be kept together, which led to both occupying the first row of the canvas. Then p3 and p4 were slotted into the second row.

Annotations

Annotations become helpful because you can add omnibus titles and tags for individual plots. For example, you can generate a common title, subtitle, caption, etc as shown below).

(p1 + p2) / (p3 + p4) +
  plot_annotation(
  title = 'The surprising truth about diamonds',
  subtitle = 'These plots will reveal untold secrets about one of our beloved data-sets',
  caption = 'Disclaimer: None of these plots are insightful',
  tag_levels = c('a', '1'),
  tag_prefix = 'Fig. ',
  tag_sep = '.',
  tag_suffix = ':'
  ) &
  theme(
    plot.tag.position = c(0, 1),
    plot.tag = element_text(size = 9, hjust = 0, vjust = 0, color = "steelblue")
    )
../_images/Module06_23_0.png

Spacing and Sizing

We can also tweak the sizes of individual rows and columns, control the space between plots, and so on. First up, spacing the plots with plot_spacer()

options(repr.plot.width = 16, repr.plot.height = 8)

(p1 + plot_spacer() + p2 + plot_spacer() + p3)
../_images/Module06_25_0.png

Sizing the plots with relative sizes?

options(repr.plot.width = 16, repr.plot.height = 16)

p1 + p2 + p3 + p4 + 
  plot_layout(widths = c(2, 1))
../_images/Module06_27_0.png

Alternatively, we could specify size with unit vectors, as shown below.

p1 + p2 + p3 + p4 + 
  plot_layout(
    widths = c(2, 1),
    heights = unit(c(5, 1), c('cm', 'null'))
    )
../_images/Module06_29_0.png

Moving Beyond the grid

We can use a layout design to get a little more flexibility but still retain full control over the result. Layout designs can be done in two ways so let us see the easiest route – as a text setup. “When using the textual representation it is your responsibility to make sure that each area is rectangular. The only exception is # which denotes empty areas and can thus be of any shape.”

layout <- "
##BBBB
AACCDD
##CCDD
"
p2 + p3 + p4 + p1 + 
  plot_layout(design = layout)
../_images/Module06_31_0.png

The other path is using area() inside layout, as shown below.

layout <- c(
  area(t = 2, l = 1, b = 5, r = 4),
  area(t = 1, l = 3, b = 3, r = 5)
  )

p3 + p4 + 
  plot_layout(design = layout)
../_images/Module06_33_0.png

Watch the specification here with wrap_plots()

layout <- '
A##
#B#
##C
'
wrap_plots(A = p1, B = p2, C = p3, design = layout)
../_images/Module06_35_0.png

Fixed-aspect plots

There are some plots that use fixed coordinates and these should not be disturbed. Here is an example where the map has fixed coordinates specified via coord_fixed(1.3)

library(urbnmapr)

ggplot() +
  geom_polygon(
    data = states, 
    aes(x = long, y = lat, group = group, fill = state_abbv)
    ) +
  coord_fixed(1.3) +
  ggthemes::theme_map() + 
  theme(legend.position = "none") +
  labs(title = "Fixed!!") -> mymap

mymap + p1 + p2 + p3
../_images/Module06_37_0.png

Mapping

Maps are very powerful visualizations because they allow you to highlight and reflect patterns, clusters, with relative ease. For example, is poverty really higher in Appalachian counties? What about the percent of the population without health insurance? Literacy? Opioid deaths; do they follow transportation routes? What about COVID-19 cases? Maps to the rescue!

Building a map requires a few elements. First and foremost, you need some data to show on a map. Second, you need to have the geographic coordinates needed to build a map, basically the latitude and longitude of the geographies (states, cities, school districts, etc.) that you want to map. Third, you want a column that contains the names of the geographies you want to map, and these should be properly formatted (i.e., in titlecase) for displaying on the map.

Let us start by building a simple state map with the {urbnmapr} package. It comes with the necessary data for states and counties, respectively, and works well with {ggplot2}. Note the reliance on geom_polygon() now.

If you get an error message about urbnmapr not found, go ahead and install it ONCE via devtools::install_github("UrbanInstitute/urbnmapr")

head(states)
A tibble: 6 × 9
longlatorderholepiecegroupstate_fipsstate_abbvstate_name
<dbl><dbl><int><lgl><fct><fct><chr><chr><chr>
-88.4732331.893861FALSE101.101ALAlabama
-88.4688831.930262FALSE101.101ALAlabama
-88.4686631.933173FALSE101.101ALAlabama
-88.4550432.039724FALSE101.101ALAlabama
-88.4549632.040585FALSE101.101ALAlabama
-88.4534232.053056FALSE101.101ALAlabama
 
head(counties)
A tibble: 6 × 12
longlatorderholepiecegroupcounty_fipsstate_abbvstate_fipscounty_namefips_classstate_name
<dbl><dbl><int><lgl><fct><fct><chr><chr><chr><chr><chr><chr>
-86.9176032.664171FALSE101001.101001AL01Autauga CountyH1Alabama
-86.8165732.660122FALSE101001.101001AL01Autauga CountyH1Alabama
-86.7133932.661733FALSE101001.101001AL01Autauga CountyH1Alabama
-86.7142232.705694FALSE101001.101001AL01Autauga CountyH1Alabama
-86.4131232.707395FALSE101001.101001AL01Autauga CountyH1Alabama
-86.4111732.409946FALSE101001.101001AL01Autauga CountyH1Alabama
 
ggplot() +
  geom_polygon(
      data = states, 
      aes(
          x = long, 
          y = lat, 
          group = group),
      fill = "white", 
      color = "steelblue"
  ) +
  coord_fixed(1.3)
../_images/Module06_42_0.png

Note the use of geom_polygon() and that longitude goes on the x-axis and latitude on the y-axis. I also specified the fill color and the border (via color = "")

Note that this is just an empty map with the shapes of the states, and also that Alaska and Hawaii have been moved so that they can be displayed on the map.

We could build a much better map by removing the x and y axis labels and tick marks, and setting a white background using theme_map() from the {ggthemes} package. We could also fill with some colors, say on the basis of the state_name.

Again, if you get an error with ggthemes, go ahead and install it ONCE via install.packages("ggthemes")

ggplot() +
  geom_polygon(
      data = states, 
      aes(
          x = long, 
          y = lat, 
          group = group, 
          fill = state_name),
      color = "white"
      ) +
  coord_fixed(1.3) +
  ggthemes::theme_map() +
  theme(legend.position = "none")
../_images/Module06_46_0.png

Note the legend has been hidden here since the legend with 50 states would take up all the space on the plotting canvas!

More importantly, this is not a very useful map because it would be much better to color the map on the basis of some substantive variable such as population density, income, crime rates, health insurance coverage rates, vaccination rates, population size, and so on. But we need data for each state to fold this information in.

Let us see what information lurks in the statedata file.

head(statedata)
A tibble: 6 × 6
yearstate_fipsstate_namehhpophoratemedhhincome
<int><chr><chr><int><dbl><int>
201501Alabama 18463800.681432944700
201502Alaska 2501830.631186070600
201504Arizona 24630120.620617851000
201505Arkansas 11446570.654603142000
201506California128954710.537221964600
201508Colorado 20745170.638942563500

Okay, two things stand out – horate (the homeownership rate), and medhhincome (the median household income). Let us fill with median household income but to do so, we will need to join statedata to our states file. Why? Because we need coordinates to map anything and statedata does not contain coordinates.

states %>%
  left_join(
      statedata, 
      by = c("state_fips", "state_name")
  ) -> state.df

left_join: added 4 columns (year, hhpop, horate, medhhincome)

           > rows only in x        0

           > rows only in y  (     0)

           > matched rows     83,933

           >                 ========

           > rows total       83,933

head(state.df)
A tibble: 6 × 13
longlatorderholepiecegroupstate_fipsstate_abbvstate_nameyearhhpophoratemedhhincome
<dbl><dbl><int><lgl><fct><fct><chr><chr><chr><int><int><dbl><int>
-88.4732331.893861FALSE101.101ALAlabama201518463800.681432944700
-88.4688831.930262FALSE101.101ALAlabama201518463800.681432944700
-88.4686631.933173FALSE101.101ALAlabama201518463800.681432944700
-88.4550432.039724FALSE101.101ALAlabama201518463800.681432944700
-88.4549632.040585FALSE101.101ALAlabama201518463800.681432944700
-88.4534232.053056FALSE101.101ALAlabama201518463800.681432944700

Now we can build the map with state.df and specify fill = medhhincome inside the aes(...) command.

options(repr.plot.width = 24, repr.plot.height = 16)

ggplot() +
  geom_polygon(
      data = state.df, 
      aes(
          x = long, 
          y = lat, 
          group = group, 
          fill = medhhincome
          ),
          color = "white"
  ) +
  coord_fixed(1.3) +
  ggthemes::theme_map() +
  labs(
    title = "Median Household Income in the States (2015)",
    fill = "Median Household Income"
        ) +
  scale_fill_viridis_c(option = "magma") +
  theme(
      legend.position = "bottom",
      legend.text = element_text(size = 14),
      legend.key.width = unit(5, 'cm'),
      title = element_text(size = 20, face = "bold")
  )
../_images/Module06_54_0.png

Notice the fill color and the corresponding legend rely upon values of medhhincome, allowing a reader to get some feel for how the low versus high median household income states cluster spatially. In this case, Maryland appears to have the highest median household income, and thn some of the New England states. West Virginia, Mississippi, and Arkansas bring up the rear.

What about working with counties instead of states? Sure, let us merge countydata with the counties file and then draw the map.

counties %>%
  left_join(
      countydata, 
      by = c("county_fips")
  ) -> county.df

left_join: added 4 columns (year, hhpop, horate, medhhincome)

           > rows only in x         0

           > rows only in y  (      0)

           > matched rows     208,874

           >                 =========

           > rows total       208,874

ggplot() +
  geom_polygon(data = county.df, 
               aes(x = long, y = lat, group = group, fill = medhhincome),
               color = "white", size = 0.05) +
  coord_fixed(1.3) +
  ggthemes::theme_map() +
  theme(
      legend.position = "bottom",
      legend.text = element_text(size = 14),
      legend.key.width = unit(5, 'cm'),
      title = element_text(size = 20, face = "bold")
  ) +
  labs(
    title = "Median Household Income in the Counties (2015)",
    fill = "Median Household Income"
        ) +
  scale_fill_viridis_c(option = "magma")
../_images/Module06_58_0.png

Maybe you are only interested in Florida?

county.df %>%
  filter(state_abbv == "FL") %>%
  ggplot() +
  geom_polygon(
               aes(x = long, y = lat, group = group,
                   fill = medhhincome),
               color = "white", size = 0.05
               ) +
  coord_fixed(1.3) +
  ggthemes::theme_map() +
  theme(
      legend.position = "bottom",
      legend.text = element_text(size = 14),
      legend.key.width = unit(5, 'cm'),
      title = element_text(size = 20, face = "bold")
  ) +
  labs(
    title = "Median Household Income in Floria Counties (2015)",
    fill = "Median Household Income"
        ) +
  scale_fill_viridis_c(option = "plasma")

filter: removed 203,791 rows (98%), 5,083 rows remaining
../_images/Module06_60_1.png

Hmm, so far so good but what if the data were for some geography not bundled with {urbnmapr}, school districts or places, for example? Not a problem, we just have to go the extra mile. First we would have to find, download, and upload the shapefile. Say I am looking for places (loosely described as municipalities) in New Hampshire. Well, the {tigris} package comes in handy because it allows you to get whatever geography’s shapefiles you want. Below I am getting the shapefile for New Hampshire.

Again, if you get an error with tigris, install it ONCE via install.packages("tigris")

library(tigris)
options(tigris_use_cache = TRUE)

places(
  state = "NH", cb = TRUE, year = 2018, progress_bar = FALSE
  ) -> places.nh

To enable 
caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.

Here, we asked for all places in New Hampsire via places(state = "NH", ...) and asked for the map boundaries to be those on file for 2018. Boundary files are updated every year or two, depending upon changes and state reporting.

places.nh %>%
    head()

Registered S3 method overwritten by 'geojsonsf':
  method        from   
  print.geojson geojson
A sf: 6 × 10
STATEFPPLACEFPPLACENSAFFGEOIDGEOIDNAMELSADALANDAWATERgeometry
<chr><chr><chr><chr><chr><chr><chr><dbl><dbl><POLYGON [°]>
13361860023780881600000US33618603361860Pittsfield 57 4114545 0POLYGON ((-71.34923 43.2990...
23367460026297381600000US33674603367460Sanbornville57 4106641 3626POLYGON ((-71.04123 43.5649...
33387140023780991600000US33871403387140Woodsville 57 2271155 63295POLYGON ((-72.04107 44.1542...
43312900008735671600000US33129003312900Claremont 251117661372323426POLYGON ((-72.41538 43.3802...
53311300023780551600000US33113003311300Charlestown 57 2107485 121062POLYGON ((-72.43447 43.2304...
63378340026297441600000US33783403378340Walpole 57 3152287 27471POLYGON ((-72.43535 43.0751...

What you have is a shapefile … a document that gives you the coordinates stored in POLYGON(...) and other information such as the unique FIPS codes, the placename, land area, water area, and so on.

Okay, so now that I have the shapefile, how can I use it?

I need to fortify it so that it looks like a regular dataframe rather than the native SpatialPolygonsDataFrame format it comes in. When I go to make the map I am going to add the state shapefile too since otherwise the state’s boundary will not show up.

places.nh %>%
  fortify(region = "GEOID") -> nh.df

names(nh.df)
  1. 'STATEFP'
  2. 'PLACEFP'
  3. 'PLACENS'
  4. 'AFFGEOID'
  5. 'GEOID'
  6. 'NAME'
  7. 'LSAD'
  8. 'ALAND'
  9. 'AWATER'
  10. 'geometry'
nh.df %>%
    head()
A sf: 6 × 10
STATEFPPLACEFPPLACENSAFFGEOIDGEOIDNAMELSADALANDAWATERgeometry
<chr><chr><chr><chr><chr><chr><chr><dbl><dbl><POLYGON [°]>
13361860023780881600000US33618603361860Pittsfield 57 4114545 0POLYGON ((-71.34923 43.2990...
23367460026297381600000US33674603367460Sanbornville57 4106641 3626POLYGON ((-71.04123 43.5649...
33387140023780991600000US33871403387140Woodsville 57 2271155 63295POLYGON ((-72.04107 44.1542...
43312900008735671600000US33129003312900Claremont 251117661372323426POLYGON ((-72.41538 43.3802...
53311300023780551600000US33113003311300Charlestown 57 2107485 121062POLYGON ((-72.43447 43.2304...
63378340026297441600000US33783403378340Walpole 57 3152287 27471POLYGON ((-72.43535 43.0751...

In the map, instead of filter(...) I am using another command to only keep boundaries for places in New Hampshire. This is being done via data = subset(state.df, state_name == "New Hampshire")

ggplot() +
    geom_polygon(
        data = subset(state.df, state_name == "New Hampshire"),
        aes(x = long, y = lat, group = group),
        fill = "white", color = "black"
        ) + 
    geom_sf(
        data = nh.df,
        aes(fill = GEOID)
        ) +
  ggthemes::theme_map() +
  theme(legend.position = "none")
../_images/Module06_71_0.png

Of course, the fill is superficial here. But say we had some data for places in New Hampshire, maybe the size of the population, as in nh.data.RData. Now we could join nh.data with nh.df to create nh and then map. Note the join keys … GEOID in each file.

load("data/nh.data.RData")

head(nh.data)
A tibble: 6 × 3
GEOIDNAMEpopulation
<chr><chr><dbl>
3300980Alton CDP, New Hampshire 168
3301220Amherst CDP, New Hampshire 709
3301620Antrim CDP, New Hampshire 1232
3301940Ashland CDP, New Hampshire 1353
3303620Bartlett CDP, New Hampshire 104
3304660Belmont CDP, New Hampshire 1814
 
nh.df %>%
  left_join(
      nh.data, 
      by = c("GEOID" = "GEOID")
  ) -> nh

left_join: added 3 columns (NAME.x, NAME.y, population)

           > rows only in x    0

           > rows only in y  ( 0)

           > matched rows     97

           >                 ====

           > rows total       97

head(nh)
A sf: 6 × 12
STATEFPPLACEFPPLACENSAFFGEOIDGEOIDNAME.xLSADALANDAWATERNAME.ypopulationgeometry
<chr><chr><chr><chr><chr><chr><chr><dbl><dbl><chr><dbl><POLYGON [°]>
13361860023780881600000US33618603361860Pittsfield 57 4114545 0Pittsfield CDP, New Hampshire 1586POLYGON ((-71.34923 43.2990...
23367460026297381600000US33674603367460Sanbornville57 4106641 3626Sanbornville CDP, New Hampshire 581POLYGON ((-71.04123 43.5649...
33387140023780991600000US33871403387140Woodsville 57 2271155 63295Woodsville CDP, New Hampshire 903POLYGON ((-72.04107 44.1542...
43312900008735671600000US33129003312900Claremont 251117661372323426Claremont city, New Hampshire 13016POLYGON ((-72.41538 43.3802...
53311300023780551600000US33113003311300Charlestown 57 2107485 121062Charlestown CDP, New Hampshire 1029POLYGON ((-72.43447 43.2304...
63378340026297441600000US33783403378340Walpole 57 3152287 27471Walpole CDP, New Hampshire 519POLYGON ((-72.43535 43.0751...

Now we plot:

ggplot() +
  geom_polygon(
    data = subset(
        state.df, 
        state_name == "New Hampshire"
    ),
    aes(
        x = long, 
        y = lat, 
        group = group
    ),
      fill = "white", 
      color = "black"
    ) +
  geom_sf(
    data = nh,
    aes(fill = population)
    ) +
  coord_sf() +
  scale_fill_viridis_c(option = "viridis") + 
  ggthemes::theme_map() +
  theme(
      legend.position = "bottom",
      legend.text = element_text(size = 14),
      legend.key.width = unit(5, 'cm'),
      title = element_text(size = 20, face = "bold")  
  ) +
  labs(
    fill = "Population Size",
    title = "Population Distribution in New Hampshire's Places",
    subtitle = "(American Community Survey, 2014-2018)"
       )
../_images/Module06_77_0.png

You could have also filled by creating quartiles, etc., using {Santoku}, so do not forget that option.

Before we move on, one more map, to show you the possibilities. Here, I am plotting the locations of the parking tickets issued in Philadelphia, this time with the {leaflet} package. First the tickets data-set, reduced to a random sample of 20% of the tickets issued in the month of December in 2017 (to keep the data size manageable). I have also created a popup that will display specific information if someone clicks on a point in the map.

Again, install any missing packages you are alerted about but remember that you need to do this ONCE, not every time you need to use the package.

readr::read_csv(
    "https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-12-03/tickets.csv"
    ) %>%
  mutate(
    year = lubridate::year(issue_datetime),
    month = lubridate::month(issue_datetime)
         ) %>%
  filter(month == 12, lon > -75.5) %>%
  sample_frac(0.2) -> tickets

Rows: 1260891 Columns: 7

── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): violation_desc, issuing_agency
dbl  (4): fine, lat, lon, zip_code
dttm (1): issue_datetime

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'month' (double) with 12 unique values and 0% NA

filter: removed 1,169,690 rows (93%), 91,201 rows remaining

sample_frac: removed 72,961 rows (80%), 18,240 rows remaining

tickets %>%
  unite(
      display, 
      c(issuing_agency, issue_datetime, fine),
      sep = "; ", 
      remove = FALSE
  ) -> tickets

library(leaflet)
library(htmltools)
library(widgetframe)

Loading required package: htmlwidgets

lst <- list()

leaflet(tickets) %>%
  addTiles() %>%
  addCircles(
      lng = ~ lon, 
      lat = ~ lat, 
      popup = ~htmlEscape(display),
      color = "steelblue", 
      opacity = 0.10
  ) -> leaf01

leaf01 -> lst

htmltools::tagList(lst)

Voila! A few lines of code and we have an interactive map that can be used to display whatever evidence we want to display. Note that you need geographic coordinates since without them the data cannot be attached to a physical location.

Let us see another twist on this. Say I am trying to map the total number of COVID-19 cases in Ohio. I know the county and I know the number of cases that occurred, as well as the latitude and longitude of each county. Well, I can build a similar plot except making the size of the circle conditional upon the number of cases. The more the cases, the larger the radius of the circle.

readr::read_csv(
    "https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv"
    ) -> covid 

Rows: 2050611 Columns: 6

── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): county, state, fips
dbl  (2): cases, deaths
date (1): date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

covid %>%
  filter(
      state == "Ohio", date == "2020-04-17"
  ) -> cov19

filter: removed 2,050,524 rows (>99%), 87 rows remaining

But there are no coordinates in cov19. Well, I can get those from the {housingData} package, and then merge that with cov19.

library(housingData)

geoCounty %>%
  filter(state == "OH") %>%
  separate(
      county,
      into = c("countyname", "extra"),
      sep = " County",
      remove = TRUE
  ) %>%
  mutate(
      countyname = stringr::str_to_sentence(countyname)
  ) -> oh

filter: removed 2,987 rows (97%), 88 rows remaining

mutate: changed one value (1%) of 'countyname' (0 new NA)

oh %>%
  left_join(
      cov19, 
      by = c("countyname" = "county")
  ) -> ohcov19

head(ohcov19)

left_join: added 7 columns (fips.x, state.x, date, state.y, fips.y, …)

           > rows only in x    2

           > rows only in y  ( 1)

           > matched rows     86

           >                 ====

           > rows total       88
A data.frame: 6 × 13
fips.xcountynameextrastate.xlonlatrMapStaterMapCountydatestate.yfips.ycasesdeaths
<fct><chr><chr><fct><dbl><dbl><fct><fct><date><chr><chr><dbl><dbl>
139001Adams OH-83.4635938.85662ohioadams 2020-04-17Ohio39001 30
239003Allen OH-84.1082540.77675ohioallen 2020-04-17Ohio39003659
339005Ashland OH-82.2692240.86122ohioashland 2020-04-17Ohio39005 50
439007AshtabulaOH-80.7564741.71017ohioashtabula2020-04-17Ohio39007544
539009Athens OH-82.0405339.34576ohioathens 2020-04-17Ohio39009 31
639011Auglaize OH-84.2219240.56421ohioauglaize 2020-04-17Ohio39011211

Now I want to create a new column that shows the county name and the number of cases.

ohcov19 %>%
  unite(
      display, 
      c(countyname, cases), 
      sep = ": ",
      remove = FALSE
  ) -> ohcov19

head(ohcov19)
A data.frame: 6 × 14
fips.xdisplaycountynameextrastate.xlonlatrMapStaterMapCountydatestate.yfips.ycasesdeaths
<fct><chr><chr><chr><fct><dbl><dbl><fct><fct><date><chr><chr><dbl><dbl>
139001Adams: 3 Adams OH-83.4635938.85662ohioadams 2020-04-17Ohio39001 30
239003Allen: 65 Allen OH-84.1082540.77675ohioallen 2020-04-17Ohio39003659
339005Ashland: 5 Ashland OH-82.2692240.86122ohioashland 2020-04-17Ohio39005 50
439007Ashtabula: 54AshtabulaOH-80.7564741.71017ohioashtabula2020-04-17Ohio39007544
539009Athens: 3 Athens OH-82.0405339.34576ohioathens 2020-04-17Ohio39009 31
639011Auglaize: 21 Auglaize OH-84.2219240.56421ohioauglaize 2020-04-17Ohio39011211

Okay, now we have everything we need to build the map.

lst <- list()

leaflet(ohcov19) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~ lon, 
    lat = ~ lat, 
    popup = ~htmlEscape(display),
    color = "salmon", 
    opacity = 0.10, 
    radius = ~sqrt(cases)
    ) -> leaf02

leaf02 -> lst

htmltools::tagList(lst)

Note that the radius of the circles is being driven by the square-root of the number of cases (a convenient transformation that makes the size managable). If we relied on the cases themselves, with no square-root, we would get this, useless map because the large number of cases would overwhelm everything else, coloring everything shades of red.

lst <- list()

leaflet(ohcov19) %>%
  addTiles() %>%
  addCircleMarkers(
      lng = ~ lon, 
      lat = ~ lat, 
      popup = ~htmlEscape(display),
      color = "salmon", 
      opacity = 0.10, 
      radius = ~cases
    ) -> leaf02bad

leaf02bad -> lst

htmltools::tagList(lst)

Interactive Graphics with Plotly and Highcharter

Interactive graphics are useful in situations where you would like the user/viewer to see the data values or other details by hovering over or clicking on the graphic. Say, for example, I have a scatterplot and want to make it interactive. How can I do that?

One crude and fast way to do that is by saving my ggplot2 object and then using {plotly} to add a ggplotly() wrapper around the plot.

In the example below I am saving the plot as pl01, then wrapping it in ggplotly with ggplotly(pl01) -> lst.

Note too that the lst <- list() and htmltools::tagList(lst) commands will show up for the interactive plots, and they need to. Otherwise the plot may not render correctly in the notebook.

library(plotly)

lst <- list()

ggplot() +
  geom_point(
    data = mpg,
    mapping = aes(
        x = cty, 
        y = hwy, 
        color = trans)
    ) +
  labs(
      x = "City Mileage",
      y = "Highway Mileage",
      color = "Transmission"
  ) -> pl01

ggplotly(pl01) -> lst

htmltools::tagList(lst)

Attaching package: ‘plotly’

The following objects are masked from ‘package:tidylog’:

    distinct, filter, group_by, mutate, rename, select, slice,
    summarise, transmute, ungroup

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

These plots are useful when presenting data to a live audience (in a talk, or on the web).

Rather than use plotly, I prefer {highcharter} since it does a lot of things well with minimal fuss, and yet the resulting plots are aesthetically pleasing.

Let us stay with the COVID-19 example. Say I want a bar-chart of the total number of cases by state and want to do this via highcharter.

library(highcharter)

covid %>%
  filter(
      date == "2020-04-17"
  ) %>%
  rename(
      State = state, 
      `Total Cases` = cases
  ) -> tab1

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

lst <- list()

hchart(
    tab1, 
    "bar", 
    hcaes(
        x = State, 
        y = `Total Cases`
        )
    ) -> hc01

hc01 -> lst
htmltools::tagList(lst)

Notice the key elements here: The basic function call is hchart() and we are specifying that we want a bar-chart, and we are also providing the quantities that should go on the x and y axis, respectively. Note that x actually ends up as the y when you specify a “bar” chart.

What if I wanted a line-chart, maybe of the number of cases over time? And I wanted this just for a few states? We could do that too, as shown below. Note that I am creating tab2, a frequency table of the number of cases by state and date, and then converting total_cases into a logarithmic form (saved as log_cases) so that we can compare the rate of change from one date to the next on a common scale.

covid %>%
  filter(
      state %in% c("Ohio", "Florida", "California", "New Jersey", "Ohio", "New York"),
      date >= "2020-03-01"
  ) %>%
  group_by(state, date) %>%
  mutate(
      log_cases = log(sum(cases))
  ) %>%
    ungroup() -> tab2

head(tab2)
A tibble: 6 × 7
datecountystatefipscasesdeathslog_cases
<date><chr><chr><chr><dbl><dbl><dbl>
2020-03-01Alameda California06001103.496508
2020-03-01Humboldt California06023103.496508
2020-03-01Los AngelesCalifornia06037103.496508
2020-03-01Marin California06041103.496508
2020-03-01Napa California06055103.496508
2020-03-01Orange California06059103.496508

There are duplicate rows per state per date because the counties remain. I will run distinct() to get rid of them.

tab2 %>%
    select(state, date, log_cases) %>%
    distinct() -> tab2_nodups

head(tab2_nodups)
A tibble: 6 × 3
statedatelog_cases
<chr><date><dbl>
California2020-03-013.4965076
Florida 2020-03-010.6931472
New York 2020-03-010.0000000
California2020-03-023.6375862
Florida 2020-03-020.6931472
New York 2020-03-020.0000000
 
lst <- list()

hchart(
    tab2_nodups, 
    "line", 
    hcaes(
        x = date, 
        y = log_cases, 
        group = state
        )
    ) -> hc02

hc02 -> lst
htmltools::tagList(lst)

Now here is a county-level chart that shows the total number of cases as of November 15, 2021.

The data are stored in tab3 created as shown below. Pay attention to this creation because we are not just creating a frequency table but also adding in a specific key we are calling code because we will need to join these data to the map data.

covid %>%
  group_by(county, state, fips) %>%
  filter(date == "2021-11-15") %>% 
  unite(
      Location, 
      c(county, state), 
      sep = ", ", 
      remove = TRUE
  ) -> tab3

head(tab3)
A grouped_df: 6 × 5
dateLocationfipscasesdeaths
<date><chr><chr><dbl><dbl>
2021-11-15Autauga, Alabama0100110407154
2021-11-15Baldwin, Alabama0100337875581
2021-11-15Barbour, Alabama01005 3648 79
2021-11-15Bibb, Alabama 01007 4317 92
2021-11-15Blount, Alabama 0100910536188
2021-11-15Bullock, Alabama01011 1523 44

The next step will be to take the counties data that has the longitude\latitude, and join it to the tab3 data. But before we can do that we will need to split the county_fips variable into two – stfips and fips. Why? Because we need to create a key that can be used to join these data to the highcharter map data. In highcharter there is a variable called code that looks like “us-al-001” for Autauga county in Alabama, and so on. Such a variable does not exist in the counties data-set.

library(urbnmapr)

data(counties)

head(counties)
A tibble: 6 × 12
longlatorderholepiecegroupcounty_fipsstate_abbvstate_fipscounty_namefips_classstate_name
<dbl><dbl><int><lgl><fct><fct><chr><chr><chr><chr><chr><chr>
-86.9176032.664171FALSE101001.101001AL01Autauga CountyH1Alabama
-86.8165732.660122FALSE101001.101001AL01Autauga CountyH1Alabama
-86.7133932.661733FALSE101001.101001AL01Autauga CountyH1Alabama
-86.7142232.705694FALSE101001.101001AL01Autauga CountyH1Alabama
-86.4131232.707395FALSE101001.101001AL01Autauga CountyH1Alabama
-86.4111732.409946FALSE101001.101001AL01Autauga CountyH1Alabama
 
counties %>%
  separate(
      county_fips, 
      into = c("stfips", "fips"), 
      sep = 2, 
      remove = FALSE
      ) %>%
  mutate(
      leader = "us", 
      stlower = stringr::str_to_lower(state_abbv)
      ) %>%
  unite(
      code, 
      c(leader, stlower, fips), 
      sep = "-"
      ) -> cdf

head(cdf)
A tibble: 6 × 14
longlatorderholepiecegroupcounty_fipsstfipscodestate_abbvstate_fipscounty_namefips_classstate_name
<dbl><dbl><int><lgl><fct><fct><chr><chr><chr><chr><chr><chr><chr><chr>
-86.9176032.664171FALSE101001.10100101us-al-001AL01Autauga CountyH1Alabama
-86.8165732.660122FALSE101001.10100101us-al-001AL01Autauga CountyH1Alabama
-86.7133932.661733FALSE101001.10100101us-al-001AL01Autauga CountyH1Alabama
-86.7142232.705694FALSE101001.10100101us-al-001AL01Autauga CountyH1Alabama
-86.4131232.707395FALSE101001.10100101us-al-001AL01Autauga CountyH1Alabama
-86.4111732.409946FALSE101001.10100101us-al-001AL01Autauga CountyH1Alabama
 
cdf %>%
  select(code, county_fips) %>%
  distinct() -> cdf2

head(cdf2)
A tibble: 6 × 2
codecounty_fips
<chr><chr>
us-al-00101001
us-al-00301003
us-al-00501005
us-al-00701007
us-al-00901009
us-al-01101011

Now we can join cdf2 to tab3 so that the code variable will become common to the resulting data-set and highcharter!

tab3 %>%
  left_join(
      cdf2, 
      by = c("fips" = "county_fips")
  ) -> tab4

left_join: added one column (code)

           > rows only in x     117

           > rows only in y  (    9)

           > matched rows     3,133

           >                 =======

           > rows total       3,250

head(tab4)
A grouped_df: 6 × 6
dateLocationfipscasesdeathscode
<date><chr><chr><dbl><dbl><chr>
2021-11-15Autauga, Alabama0100110407154us-al-001
2021-11-15Baldwin, Alabama0100337875581us-al-003
2021-11-15Barbour, Alabama01005 3648 79us-al-005
2021-11-15Bibb, Alabama 01007 4317 92us-al-007
2021-11-15Blount, Alabama 0100910536188us-al-009
2021-11-15Bullock, Alabama01011 1523 44us-al-011

Here comes the map!

Note that we are asking forthe cases column to be used for the values that will color each county, and we are asking the map file be joined with hc-key in the highcharter file and code in tab4. The county birders will be in steelblue, and there will be 10 values used to create the fill color palette. The legend will be aligned right, and horizontal.

library(viridis)

lst <- list()

hcmap("countries/us/us-all-all", 
      data = tab4,
      name = "COVID-19 Cases", value = "cases",
      joinBy = c("hc-key", "code"),
      borderColor = "steelblue") %>%
  hc_colorAxis(stops = color_stops(10, rev(magma(10)))) %>% 
  hc_legend(layout = "horizontal", align = "right",
            floating = TRUE, valueDecimals = 0, valueSuffix = ""
           ) -> hc03

hc03 -> lst

htmltools::tagList(lst)

Loading required package: viridisLite

Note that countries/us/us-all-all indicates that we want counties. If we wanted the states instead it would have been countries/us/us-all.

What if we wanted only Ohio?

Well, in that case we could subset as shown below. In particular, we are asking that if R sees the string “oh” in a variable called code, it should keep only these rows, and save these filtered rows of data as tab5.

tab4 %>%
  filter(
      grepl("-oh-", code)
  ) -> tab5

head(tab5)
A grouped_df: 6 × 6
dateLocationfipscasesdeathscode
<date><chr><chr><dbl><dbl><chr>
2021-11-15Adams, Ohio 39001 4360105us-oh-001
2021-11-15Allen, Ohio 3900317443312us-oh-003
2021-11-15Ashland, Ohio 39005 7413146us-oh-005
2021-11-15Ashtabula, Ohio3900711166217us-oh-007
2021-11-15Athens, Ohio 39009 8083 93us-oh-009
2021-11-15Auglaize, Ohio 39011 7334110us-oh-011
 
lst <- list()

hcmap("countries/us/us-oh-all", 
      data = tab5,
      name = "COVID-19 Cases", value = "cases",
      joinBy = c("hc-key", "code"),
      borderColor = "steelblue") %>%
  hc_colorAxis(stops = color_stops(10, rev(magma(10)))) %>% 
  hc_legend(layout = "horizontal", align = "right",
            floating = TRUE, valueDecimals = 0, valueSuffix = ""
           ) -> hc04

hc04 -> lst

htmltools::tagList(lst)

There you have it!

The one downside to these interactive charts is that they are best displayed in html files but in PDF and Word document they lose that interactivity. Hence you see them a lot on blogs and other web-based documents.

All of these packages have been growing so it is quite likely that as software development continues even that barrier might be eliminated.

Exercises for Practice

Exercise 01

Create a map of all the counties in New York. Be sure to title the map and to fill in each county with the total number of COVID19 cases they have seen to date. In addition, draw county borders in white. Use theme_map() and make sure the legend is at the bottom. [Hint: You will need to calculate the total number of cases per county and then join the resulting file with the counties data file to get the latitude/longitudes for the counties.]

Exercise 02

Run the following code chunk to load data on the murder, assault and rape rates per 100,000 persons. Urbanpop is the percent of the state population that lives in an urban area.

library(tidyverse)

data(USArrests)
names(USArrests)
USArrests$statename <- rownames(USArrests)

head(USArrests)
  1. 'Murder'
  2. 'Assault'
  3. 'UrbanPop'
  4. 'Rape'
A data.frame: 6 × 5
MurderAssaultUrbanPopRapestatename
<dbl><int><int><dbl><chr>
Alabama13.22365821.2Alabama
Alaska10.02634844.5Alaska
Arizona 8.12948031.0Arizona
Arkansas 8.81905019.5Arkansas
California 9.02769140.6California
Colorado 7.92047838.7Colorado

Now create a state-level map of the 50 states making sure to use UrbanPop to fill each state. Title the map and place the legend at the bottom.

Exercise 03

Use the USArrests data to draw scatterplots of (a) Murder versus UrbanPop, (b) Assault versus UrbanPop, and (c) Rape versus UrbanPop. Save each of these scatterplots by name and then use patchwork to create a single canvas that includes all three plots. Make sure you label the x-axis, y-axis, and title each plot.

Exercise 04

Now create highcharter versions of each of the three scatterplots you created in Exercise (3) above. You should end up with three scatterplots, each on its own canvas.